Mathematica for Quantum Information Theory

The basics

This notebook was inspired by Jony Hudson’s Mathematica course (which is now rather hard to find online!)

We will quickly revisit the Pauli matrices from week 1:

Section 1: Pauli matrices and single qubit gates

Let’s define the Pauli matrices as

Id ={{1, 0}, {0, 1}};
X = {{0, 1}, {1, 0}};
Y = {{0, -I}, {I, 0}};
Z = {{1, 0}, {0, -1}};

paulis = {X, Y, Z};

We naughtily defined these as upper-case letters, which you should generally avoid doing so as not to collide with inbuilt functions like N. We could have done this much more concisely using the PauliMatrix function:

paulis = {X,Y,Z} = PauliMatrix[Range[3]];

To select the kth entry of a list X in Mathematica we can use X[[k]]. Note the indices start at 1 (not 0)! To print matrices nicely in Mathematica we can then use the MatrixForm function. e.g.

paulis[[2]] // MatrixForm

where // is just the postfix syntax for applying a function:

f @ a

a // f

Note that directly multiplying matrices...

X Y

X * Y

will multiply them element-wise, i.e. where (Subscript[(XY), ij]) = (Subscript[X, ii])(Subscript[Y, ij]). Matrix multiplication as we ordinarily define it is performed using Dot:

X.Y //MatrixForm

Let’s define a function to compute the commutation relation between two matrices.

commutator[A_, B_]:= A.B - B.A

commutator[X,Y] // MatrixForm

The operator := is SetDelayed, which you should typically always use (instead of just =) when defining functions. The latter evaluations the function definition immediately which often isn’t what we want!

Define a function that computes the anti-commutator and produce a Table (or take a more sophisticated approach if you like) to verify that the standard anti-commutation rules hold

Solution

anticommutator[A_, B_]:= A.B + B.A

All anti-commutators are:

Table[
	anticommutator[a,b] // MatrixForm,
	{a,paulis}, {b,paulis}
]

We can make this prettier and clearer by using StringJoin, Row and Column:

"a" <> "b"

Row[{1, A, "hello", π/.5}]

Column[{1, A, (({{a, b}, {c, d}})), Plot[x,{x,0,1}]}]

to produce

label[X] = "X";
label[Y] = "Y";
label[Z] = "Z";
label[a_,b_] := "{" <> label[a] <> "," <> label[b] <> "} =";

table = Table[
	Row[{label[a,b],anticommutator[a,b] // MatrixForm}],
	{a,paulis}, {b,paulis}
];

Column @ Flatten @ table

Define a function that implements evolution under a Pauli operator by some angle. That is, create a function that evaluates (e^(-i θ P)) where θ is a scalar and P is either X, Y or Z

Hints:
1) Use MatrixExp.
2) Mathematica has a long list of helpful short cuts for symbols
e.g. see https://reference.wolfram.com/language/guide/MathematicalTypesetting.html and https://reference.wolfram.com/language/guide/SpecialCharacters.html

Solution

rotation[theta_, sigma_]:= MatrixExp[-I theta sigma]

rotation[θ, X] // MatrixForm

rotation[θ, Y] // MatrixForm

rotation[θ, Z] // MatrixForm

Generalize this to implement evolution about an arbitrary axis r on the Bloch Sphere, i.e. (e^(-i θ n .σ)) where n is any 3D vector and σ = {X, Y, Z}

Solution

rotation3D[theta_, n_]:= MatrixExp[-I theta  n.paulis]

rotation3D[θ, {1,0,0}] // MatrixForm
rotation3D[θ, {0,1,0}] // MatrixForm
rotation3D[θ, {0,0,1}] // MatrixForm

rotation3D[θ, {1,2,3}] // MatrixForm

rotation3D[θ, {nx,ny,nz}] // MatrixForm

Simplify the expression for (e^(-i θ n .σ)) when n is a real unit-vector.

Solution

FullSimplify[
	rotation3D[θ, {nx,ny,nz}],
	(nx^2)+(ny^2)+(nz^2)==1
] // MatrixForm

Define a function that implements the standard analytic expression Pauli exponentiation, i.e.  (e^(-iθ (OverHat[n]).σ)) = cos(θ) I - i sin(θ)  (OverHat[n]) . σ and use this to check that your function above is correct.

Solution

euler3D[theta_, n_] := Cos[theta] IdentityMatrix[2] - I Sin[theta] n . paulis

If these expressions are equivalent, their difference must be the zero matrix:

euler3D[θ, {nx, ny,nz}]

FullSimplify[
	rotation3D[θ,{nx,ny,nz}] - euler3D[θ,{nx,ny,nz}],
	(nx^2)+(ny^2)+(nz^2)==1
]

Assessed problem

Define a function to create an arbitrary single qubit state and hence convince yourself that evolution under (e^(i α (OverHat[n]).σ)) implements a rotation by an angle 2α about unit vector (OverHat[n]) on the Bloch sphere.

Hint: Use the representation of the Bloch vector in spherical coordinates

Solution

(* https://www.youtube.com/watch?v=7b4ULg131PM *)

Section 2: Multi-qubit systems and gates

Define projectors onto the 0 and 1 states of a qubit. Use these, and the KroneckerProduct function, to define the CNOT gate (controlling upon the left qubit)

Solution

proj0 = {{1, 0}, {0, 0}};
proj1 = {{0, 0}, {0, 1}};

CNOT = KroneckerProduct[proj0, Id] + KroneckerProduct[proj1,X]; 
CNOT // MatrixForm

The Hadamard gate is defined as follows

H = (1/(Sqrt[2])){{1,1},{1,-1}};
H // MatrixForm

It is used to prepare the plus and minus states, e.g.

H.{1,0}
H.{0,1}

Show that the Hadamard gates and CNOT gates can be combined to prepare a Bell state. (How would you prepare the other Bell states?)

Solution

We use the Hadamard to create superposition of the qubit which is subsequently controlled-upon in the CNOT gate (the left qubit). Operators act via multiplication to the right, so the CNOT appears on the left:

makeBell = CNOT . KroneckerProduct[H,Id];

The unitary acts upon the zero state |00⟩ = {1,0,0,0} to produce the Bell state (1/(Sqrt[2]))(|00⟩+|11⟩):

makeBell . {1,0,0,0}

The other Bell states are prepared by applying the same unitary upon the other possible classical input states, i.e. |01⟩, |10⟩ and |11⟩

makeBell . {0,1,0,0}
makeBell . {0,0,1,0}
makeBell . {0,0,0,1}

Below, we define calcPartialTrace which accepts a matrix and a list of qubits (indexed from zero).

getPartialKet[bitInds_List, bitVals_List, numAllBits_] := Module[
	{ops = ConstantArray[IdentityMatrix[2], numAllBits]},
	Do[
		ops[[ bitInds[[i]]+1 ]] = If[bitVals[[i]]===0, {{1},{0}}, {{0},{1}}],
		{i, Length@bitInds}];
	KroneckerProduct @@ Reverse[ops]]

calcPartialTrace[ρ_?MatrixQ, qubits_List] /; Length@qubits < Log2@Length@ρ :=
	Sum[
		With[
			{ket = getPartialKet[qubits, vals, Log2@Length@ρ]},
			ConjugateTranspose[ket] . ρ . ket],
		{vals, Tuples[{0,1}, Length@qubits]}]

It can be used like so:

ρm = Table[ρ[i,j],{i,4},{j,4}];
ρm // MatrixForm

calcPartialTrace[ρm, {0}] // MatrixForm

χ = Table[RandomComplex[],(2^3),(2^3)];
ρm = χ . (χ^\[ConjugateTranspose])/Tr[χ . (χ^\[ConjugateTranspose])] // Chop;
ρm // MatrixForm

calcPartialTrace[ρm, {0,2}] // Chop // MatrixForm

Show that tracing out half of any Bell state generates the maximally mixed state

Solution

We must first produce a density matrix from a given bell state, via ρ = |ψ⟩⟨ψ|, and then pass it to calcPartialTrace

getDensityMatrix[state_] := With[
	{colvec = List /@ state},
	colvec . ConjugateTranspose[colvec]
]

traceHalf[state_, qubit_] := 
	calcPartialTrace[getDensityMatrix[state], {qubit}]

maximallyMixed = (1/2)IdentityMatrix[2]

As expected, tracing out qubit 0 of the Bell state produces the maximally mixed state:

traceHalf[makeBell . {1,0,0,0}, 0]

We confirm it for all Bell states and all qubits:

Table[
	in = UnitVector[4,instate];
	out = traceHalf[makeBell . in, qubit];
	out == maximallyMixed,
	{instate, 1,4}, {qubit, 0,1}
]

Assessed problem

Suppose you generate a generic entangled state with the following circuit:



Write a function to calculate the reduced state of each qubit given an arbitrary angle θ and unit vector n. Does it matter which qubit you look at?

Solution

(* not telling *)

Use Manipulate and Plot to study how the purity of the reduced state depends on the rotation angle ϑ  and the rotation direction r. Comment on what you observe. 
(Hint: think about the best way to represent r)

Solution

(* c: *)